n_lambda = 10;
lambda_grid = linspace(0.01,0.15,n_lambda);
A_star_grid = nan(n_lambda,T-1);
S_star_grid = nan(n_lambda,T-1);

for i = 1:n_lambda
    disp(i)
    lambda = lambda_grid(i);
    A_star = find_A_star(tau_q,tau_e,alpha,T,lambda,phi,y_high,y_low,p);
    A_star_grid(i,:) = A_star;
    for t = 1:T-1
        A = A_star(t);
        S = surplus(tau_q,tau_e,A,alpha);
        S_star_grid(i,t) = S(2);
    end 
end